Class 23

23 Bayesian data analysis 3

23.1 Load packages and set plotting theme

library("knitr")      # for knitting RMarkdown 
library("kableExtra") # for making nice tables
library("janitor")    # for cleaning column names
library("tidybayes")  # tidying up results from Bayesian models
library("brms")       # Bayesian regression models with Stan
library("patchwork")  # for making figure panels
library("gganimate")  # for animations
library("GGally")     # for pairs plot
library("bayesplot")  # for visualization of Bayesian model fits 
library("tidyverse")  # for wrangling, plotting, etc. 

theme_set(theme_classic() + # set the theme 
            theme(text = element_text(size = 20))) # set the default text size

23.2 Load data set

Load the poker data set.

df.poker = read_csv("data/poker.csv") %>% 
  mutate(skill = factor(skill,
                        levels = 1:2,
                        labels = c("expert", "average")),
         skill = fct_relevel(skill, "average", "expert"),
         hand = factor(hand,
                       levels = 1:3,
                       labels = c("bad", "neutral", "good")),
         limit = factor(limit,
                        levels = 1:2,
                        labels = c("fixed", "none")),
         participant = 1:n()) %>% 
  select(participant, everything())

23.3 Poker

23.3.1 Visualization

Let’s visualize the data first:

df.poker %>% 
  ggplot(mapping = aes(x = hand,
                       y = balance,
                       fill = hand)) + 
  geom_point(alpha = 0.2,
             position = position_jitter(height = 0, width = 0.1)) + 
  stat_summary(fun.data = "mean_cl_boot",
               geom = "linerange",
               size = 1) + 
  stat_summary(fun.y = "mean",
               geom = "point",
               shape = 21,
               size = 4) +
  labs(y = "final balance (in Euros)") + 
  scale_fill_manual(values = c("red", "orange", "green")) +
  theme(legend.position = "none")

23.3.2 Linear model

And let’s now fit a simple (frequentist) regression model:

fit.lm = lm(formula = balance ~ 1 + hand,
            data = df.poker)

fit.lm %>% summary()

#> 
#> Call:
#> lm(formula = balance ~ 1 + hand, data = df.poker)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -12.9264  -2.5902  -0.0115   2.6573  15.2834 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   5.9415     0.4111  14.451  < 2e-16 ***
#> handneutral   4.4051     0.5815   7.576 4.55e-13 ***
#> handgood      7.0849     0.5815  12.185  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 4.111 on 297 degrees of freedom
#> Multiple R-squared:  0.3377, Adjusted R-squared:  0.3332 
#> F-statistic:  75.7 on 2 and 297 DF,  p-value: < 2.2e-16

23.3.3 Bayesian model

Now, let’s fit a Bayesian regression model using the brm() function:

fit.brm1 = brm(formula = balance ~ 1 + hand,
               data = df.poker,
               file = "cache/brm1")

fit.brm1 %>% summary()

#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: balance ~ 1 + hand 
#>    Data: df.poker (Number of observations: 300) 
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 4000
#> 
#> Population-Level Effects: 
#>             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept       5.93      0.41     5.12     6.72 1.00     2986     2744
#> handneutral     4.41      0.58     3.30     5.55 1.00     3497     2903
#> handgood        7.10      0.58     5.99     8.29 1.00     3545     2932
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     4.12      0.17     3.81     4.46 1.00     3650     2921
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
#> is a crude measure of effective sample size, and Rhat is the potential 
#> scale reduction factor on split chains (at convergence, Rhat = 1).

I use the file = argument to save the model’s results so that when I run this code chunk again, the model doesn’t need to be fit again (fitting Bayesian models takes a while …).

23.3.3.1 Visualize the posteriors

Let’s visualize what the posterior for the different parameters looks like. We use the geom_halfeyeh() function from the “tidybayes” package to do so:

fit.brm1 %>% 
  posterior_samples() %>% 
  select(-lp__) %>% 
  gather("variable", "value") %>% 
  ggplot(data = .,
         mapping = aes(y = variable, x = value)) +
  geom_halfeyeh()

And let’s look at how the samples from the posterior are correlated with each other:

fit.brm1 %>% 
  posterior_samples() %>% 
  select(b_Intercept:sigma) %>% 
  ggpairs(lower = list(continuous = wrap("points", alpha = 0.03)),
          upper = list(continuous = wrap("cor", size = 6))) + 
  theme(panel.grid.major = element_blank(),
        text = element_text(size = 12))

23.3.3.2 Compute highest density intervals

To compute the MAP (maximum a posteriori probability) estimate and highest density interval, we use the mode_hdi() function that comes with the “tidybayes” package.

fit.brm1 %>% 
  posterior_samples() %>% 
  clean_names() %>% 
  select(starts_with("b_"), sigma) %>% 
  mode_hdi() %>% 
  pivot_longer(cols = -c(.width:.interval),
               names_to = "index",
               values_to = "value") %>% 
  select(index, value) %>% 
  mutate(index = ifelse(str_detect(index, fixed(".")), index, str_c(index, ".mode"))) %>% 
  separate(index, into = c("parameter", "type"), sep = "\\.") %>% 
  pivot_wider(names_from = type, 
              values_from = value) %>% 
  kable(digits = 2) %>% 
  kable_styling(bootstrap_options = "striped",
                full_width = F)

parameter	mode	lower	upper
b_intercept	6.02	5.13	6.73
b_handneutral	4.46	3.29	5.54
b_handgood	7.07	5.95	8.21
sigma	4.11	3.81	4.46

23.3.3.3 Posterior predictive check

To check whether the model did a good job capturing the data, we can simulate what future data the Baysian model predicts, now that it has learned from the data we feed into it.

pp_check(fit.brm1, nsamples = 100)

This looks good! The predicted shaped of the data based on samples from the posterior distribution looks very similar to the shape of the actual data.

Let’s make a hypothetical outcome plot that shows what concrete data sets the model would predict:

# generate predictive samples 
df.predictive_samples = fit.brm1 %>% 
  posterior_samples() %>% 
  clean_names() %>% 
  select(contains("b_"), sigma) %>% 
  sample_n(size = 20) %>% 
  mutate(sample = 1:n()) %>% 
  group_by(sample) %>% 
  nest() %>% 
  mutate(bad = map(data, ~ .$b_intercept + rnorm(100, sd = .$sigma)),
         neutral = map(data, ~ .$b_intercept + .$b_handneutral + rnorm(100, sd = .$sigma)),
         good = map(data, ~ .$b_intercept + .$b_handgood + rnorm(100, sd = .$sigma))) %>% 
  unnest(c(bad, neutral, good)) %>% 
  select(-data)

# plot the results as an animation
p = df.predictive_samples %>% 
  pivot_longer(cols = -sample,
               names_to = "hand",
               values_to = "balance") %>% 
  mutate(hand = factor(hand, levels = c("bad", "neutral", "good"))) %>% 
  ggplot(mapping = aes(x = hand,
                       y = balance,
                       fill = hand)) + 
  geom_point(alpha = 0.2,
             position = position_jitter(height = 0, width = 0.1)) + 
  stat_summary(fun.data = "mean_cl_boot",
               geom = "linerange",
               size = 1) + 
  stat_summary(fun.y = "mean",
               geom = "point",
               shape = 21,
               size = 4) +
  labs(y = "final balance (in Euros)") + 
  scale_fill_manual(values = c("red", "orange", "green")) +
  theme(legend.position = "none") + 
  transition_manual(sample)

animate(p, nframes = 120, width = 800, height = 600, res = 96, type = "cairo")

# anim_save("poker_posterior_predictive.gif")

23.3.3.4 Test hypothesis

One key advantage of Bayesian over frequentist analysis is that we can test hypothesis in a very flexible manner by directly probing our posterior samples in different ways.

We may ask, for example, what the probability is that the parameter for the difference between a bad hand and a neutral hand (b_handneutral) is greater than 0. Let’s plot the posterior distribution together with the criterion:

fit.brm1 %>% 
  posterior_samples() %>% 
  select(b_handneutral) %>% 
  pivot_longer(cols = everything(),
               names_to = "variable",
               values_to = "value") %>% 
  ggplot(data = .,
         mapping = aes(y = variable, x = value)) +
  geom_halfeyeh() + 
  geom_vline(xintercept = 0,
             color = "red")

We see that the posterior is definitely greater than 0.

We can ask many different kinds of questions about the data by doing basic arithmetic on our posterior samples. The hypothesis() function makes this even easier. Here are some examples:

# the probability that the posterior for handneutral is less than 0
hypothesis(fit.brm1,
           hypothesis = "handneutral < 0")

#> Hypothesis Tests for class b:
#>          Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
#> 1 (handneutral) < 0     4.41      0.58     3.48     5.35          0         0
#>   Star
#> 1     
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

# the probability that the posterior for handneutral is greater than 4
hypothesis(fit.brm1,
           hypothesis = "handneutral > 4")

#> Hypothesis Tests for class b:
#>              Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1 (handneutral)-(4) > 0     0.41      0.58    -0.52     1.35       3.05
#>   Post.Prob Star
#> 1      0.75     
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

# the probability that good hands make twice as much as bad hands
hypothesis(fit.brm1,
           hypothesis = "Intercept + handgood > 2 * Intercept")

#> Hypothesis Tests for class b:
#>                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1 (Intercept+handgo... > 0     1.18      0.93    -0.29     2.67        8.5
#>   Post.Prob Star
#> 1      0.89     
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

# the probability that neutral hands make less than the average of bad and good hands
hypothesis(fit.brm1,
           hypothesis = "Intercept + handneutral < (Intercept + Intercept + handgood) / 2")

#> Hypothesis Tests for class b:
#>                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1 (Intercept+handne... < 0     0.86       0.5     0.05     1.68       0.04
#>   Post.Prob Star
#> 1      0.04     
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

Let’s double check one example, and calculate the result directly based on the posterior samples:

df.hypothesis = fit.brm1 %>% 
  posterior_samples() %>% 
  clean_names() %>% 
  select(starts_with("b_")) %>% 
  mutate(neutral = b_intercept + b_handneutral,
         bad_good_average = (b_intercept + b_intercept + b_handgood)/2,
         hypothesis = neutral < bad_good_average)

df.hypothesis %>% 
  summarize(p = sum(hypothesis)/n())

#>         p
#> 1 0.03975

23.3.3.5 Bayes factor

Another way of testing hypothesis is via the Bayes factor. Let’s fit the two models we are interested in comparing with each other:

fit.brm2 = brm(formula = balance ~ 1 + hand,
               data = df.poker,
               save_all_pars = T,
               file = "cache/brm2")

fit.brm3 = brm(formula = balance ~ 1 + hand + skill,
               data = df.poker,
               save_all_pars = T,
               file = "cache/brm3")

And then compare the models useing the bayes_factor() function:

bayes_factor(fit.brm3, fit.brm2)

#> Iteration: 1
#> Iteration: 2
#> Iteration: 3
#> Iteration: 4
#> Iteration: 5
#> Iteration: 1
#> Iteration: 2
#> Iteration: 3
#> Iteration: 4

#> Estimated Bayes factor in favor of bridge1 over bridge2: 3.83284

23.3.3.6 Full specification

So far, we have used the defaults that brm() comes with and not bothered about specifiying the priors, etc.

23.3.3.6.1 Getting the priors

Notice that we didn’t specify any priors in the model. By default, “brms” assigns weakly informative priors to the parameters in the model. We can see what these are by running the following command:

fit.brm1 %>% 
  prior_summary()

#>                  prior     class        coef group resp dpar nlpar bound
#> 1                              b                                        
#> 2                              b    handgood                            
#> 3                              b handneutral                            
#> 4 student_t(3, 10, 10) Intercept                                        
#> 5  student_t(3, 0, 10)     sigma

We can also get information about which priors need to be specified before fitting a model:

get_prior(formula = balance ~ 1 + hand,
          family = "gaussian",
          data = df.poker)

#>                  prior     class        coef group resp dpar nlpar bound
#> 1                              b                                        
#> 2                              b    handgood                            
#> 3                              b handneutral                            
#> 4 student_t(3, 10, 10) Intercept                                        
#> 5  student_t(3, 0, 10)     sigma

Here is an example for what a more complete model specification could look like:

fit.brm4 = brm(formula = balance ~ 1 + hand,
               family = "gaussian",
               data = df.poker,
               prior = c(prior(normal(0, 10), class = "b", coef = "handgood"),
                         prior(normal(0, 10), class = "b", coef = "handneutral"),
                         prior(student_t(3, 3, 10), class = "Intercept"),
                         prior(student_t(3, 0, 10), class = "sigma")),
               inits = list(list(Intercept = 0, sigma = 1, handgood = 5, handneutral = 5),
                            list(Intercept = -5, sigma = 3, handgood = 2, handneutral = 2),
                            list(Intercept = 2, sigma = 1, handgood = -1, handneutral = 1),
                            list(Intercept = 1, sigma = 2, handgood = 2, handneutral = -2)),
               iter = 4000,
               warmup = 1000,
               chains = 4,
               file = "cache/brm4",
               seed = 1)

fit.brm4 %>% summary()

#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: balance ~ 1 + hand 
#>    Data: df.poker (Number of observations: 300) 
#> Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 12000
#> 
#> Population-Level Effects: 
#>             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept       5.96      0.41     5.15     6.76 1.00     9139     8795
#> handneutral     4.37      0.58     3.23     5.53 1.00     9638     8040
#> handgood        7.05      0.58     5.93     8.19 1.00     9782     8367
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     4.13      0.17     3.81     4.49 1.00    12982     9331
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
#> is a crude measure of effective sample size, and Rhat is the potential 
#> scale reduction factor on split chains (at convergence, Rhat = 1).

We can also take a look at the Stan code that the brm() function creates:

fit.brm4 %>% stancode()

#> // generated with brms 2.7.0
#> functions { 
#> } 
#> data { 
#>   int<lower=1> N;  // total number of observations 
#>   vector[N] Y;  // response variable 
#>   int<lower=1> K;  // number of population-level effects 
#>   matrix[N, K] X;  // population-level design matrix 
#>   int prior_only;  // should the likelihood be ignored? 
#> } 
#> transformed data { 
#>   int Kc = K - 1; 
#>   matrix[N, K - 1] Xc;  // centered version of X 
#>   vector[K - 1] means_X;  // column means of X before centering 
#>   for (i in 2:K) { 
#>     means_X[i - 1] = mean(X[, i]); 
#>     Xc[, i - 1] = X[, i] - means_X[i - 1]; 
#>   } 
#> } 
#> parameters { 
#>   vector[Kc] b;  // population-level effects 
#>   real temp_Intercept;  // temporary intercept 
#>   real<lower=0> sigma;  // residual SD 
#> } 
#> transformed parameters { 
#> } 
#> model { 
#>   vector[N] mu = temp_Intercept + Xc * b;
#>   // priors including all constants 
#>   target += normal_lpdf(b[1] | 0, 10); 
#>   target += normal_lpdf(b[2] | 0, 10); 
#>   target += student_t_lpdf(temp_Intercept | 3, 3, 10); 
#>   target += student_t_lpdf(sigma | 3, 0, 10)
#>     - 1 * student_t_lccdf(0 | 3, 0, 10); 
#>   // likelihood including all constants 
#>   if (!prior_only) { 
#>     target += normal_lpdf(Y | mu, sigma);
#>   } 
#> } 
#> generated quantities { 
#>   // actual population-level intercept 
#>   real b_Intercept = temp_Intercept - dot_product(means_X, b); 
#> }

One thing worth noticing: by default, “brms” centers the predictors which makes it easier to assign a default prior over the intercept.

23.3.3.7 Inference diagnostics

So far, we’ve assumed that the inference has worked out. We can check this by running plot() on our brm object:

plot(fit.brm1)

Let’s make our own version of a trace plot for one parameter in the model:

fit.brm1 %>% 
  spread_draws(b_Intercept) %>% 
  clean_names() %>% 
  mutate(chain = as.factor(chain)) %>% 
  ggplot(aes(x = iteration, y = b_intercept, group = chain, color = chain)) + 
  geom_line()

We can also take a look at the auto-correlation plot. Ideally, we want to generate independent samples from the posterior. So we don’t want subsequent samples to be strongly correlated with each other. Let’s take a look:

variables = fit.brm1 %>%
  get_variables() %>%
  .[1:4]

fit.brm1 %>% 
  posterior_samples() %>% 
  mcmc_acf(pars = variables,
           lags = 4)

Looking good! The autocorrelation should become very small as the lag increases (indicating that we are getting independent samples from the posterior).

23.3.3.7.1 When things go wrong

Let’s try to fit a model to very little data (just two observations) with extremely uninformative priors:

df.data = tibble(y = c(-1, 1))

fit.brm5 = brm(data = df.data,
               family = gaussian,
               formula = y ~ 1,
               prior = c(prior(uniform(-1e10, 1e10), class = Intercept),
                         prior(uniform(0, 1e10), class = sigma)),
               inits = list(list(Intercept = 0, sigma = 1),
                            list(Intercept = 0, sigma = 1)),
               iter = 4000,
               warmup = 1000,
               chains = 2,
               file = "cache/brm5")

Let’s take a look at the posterior distributions of the model parameters:

summary(fit.brm5)

#> Warning: The model has not converged (some Rhats are > 1.1). Do not analyse the results! 
#> We recommend running more iterations and/or setting stronger priors.

#> Warning: There were 1203 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
#> See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: y ~ 1 
#>    Data: df.data (Number of observations: 2) 
#> Samples: 2 chains, each with iter = 4000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 6000
#> 
#> Population-Level Effects: 
#>               Estimate     Est.Error       l-95% CI      u-95% CI Rhat Bulk_ESS
#> Intercept 357550121.58 1416057299.71 -2244033111.47 3333594132.43 1.78        3
#>           Tail_ESS
#> Intercept       24
#> 
#> Family Specific Parameters: 
#>            Estimate     Est.Error l-95% CI      u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 1524412740.64 2392424321.98 21668.93 8317582240.06 1.40        4       41
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
#> is a crude measure of effective sample size, and Rhat is the potential 
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Not looking good – The estimates and credible intervals are off the charts. And the effective samples sizes in the chains are very small.

Let’s visualize the trace plots:

plot(fit.brm5)

fit.brm5 %>% 
  spread_draws(b_Intercept) %>% 
  clean_names() %>% 
  mutate(chain = as.factor(chain)) %>% 
  ggplot(aes(x = iteration,
             y = b_intercept,
             group = chain,
             color = chain)) + 
  geom_line()

Given that we have so little data in this case, we need to help the model a little bit by providing some slighlty more specific priors.

fit.brm6 = brm(data = df.data,
               family = gaussian,
               formula = y ~ 1,
               prior = c(prior(normal(0, 10), class = Intercept), # more reasonable priors
                         prior(cauchy(0, 1), class = sigma)),
               iter = 4000,
               warmup = 1000,
               chains = 2,
               seed = 1,
               file = "cache/brm6")

Let’s take a look at the posterior distributions of the model parameters:

summary(fit.brm6)

#> Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
#> See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: y ~ 1 
#>    Data: df.data (Number of observations: 2) 
#> Samples: 2 chains, each with iter = 4000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 6000
#> 
#> Population-Level Effects: 
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept    -0.13      1.69    -4.10     3.06 1.00     1352      909
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     2.04      1.88     0.61     6.95 1.00     1144     1426
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
#> is a crude measure of effective sample size, and Rhat is the potential 
#> scale reduction factor on split chains (at convergence, Rhat = 1).

This looks much better. There is still quite a bit of uncertainty in our paremeter estimates, but it has reduced dramatically.

Let’s visualize the trace plots:

plot(fit.brm6)

fit.brm6 %>% 
  spread_draws(b_Intercept, sigma) %>% 
  clean_names() %>% 
  mutate(chain = as.factor(chain)) %>% 
  pivot_longer(cols = c(b_intercept, sigma)) %>% 
  ggplot(aes(x = iteration,
             y = value,
             group = chain,
             color = chain)) + 
  geom_line() + 
  facet_wrap(vars(name), ncol = 1)

Looking mostly good!

23.4 Additional resources

• Tutorial on visualizing brms posteriors with tidybayes
• Hypothetical outcome plots
• Visual MCMC diagnostics
• How to model slider data the Baysian way

23.5 Session info

Information about this R session including which version of R was used, and what packages were loaded.

sessionInfo()

#> R version 3.6.2 (2019-12-12)
#> Platform: x86_64-apple-darwin15.6.0 (64-bit)
#> Running under: macOS Mojave 10.14.6
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] forcats_0.4.0    stringr_1.4.0    dplyr_0.8.3      purrr_0.3.3     
#>  [5] readr_1.3.1      tidyr_1.0.0      tibble_2.1.3     tidyverse_1.3.0 
#>  [9] bayesplot_1.7.1  GGally_1.4.0     gganimate_1.0.4  ggplot2_3.2.1   
#> [13] patchwork_1.0.0  brms_2.10.0      Rcpp_1.0.3       tidybayes_1.1.0 
#> [17] janitor_1.2.0    kableExtra_1.1.0 knitr_1.26      
#> 
#> loaded via a namespace (and not attached):
#>   [1] readxl_1.3.1              backports_1.1.5          
#>   [3] Hmisc_4.3-0               plyr_1.8.5               
#>   [5] igraph_1.2.4.2            lazyeval_0.2.2           
#>   [7] splines_3.6.2             svUnit_0.7-12            
#>   [9] crosstalk_1.0.0           rstantools_2.0.0         
#>  [11] inline_0.3.15             digest_0.6.25            
#>  [13] htmltools_0.4.0           rsconnect_0.8.16         
#>  [15] fansi_0.4.0               checkmate_1.9.4          
#>  [17] magrittr_1.5              cluster_2.1.0            
#>  [19] modelr_0.1.5              matrixStats_0.55.0       
#>  [21] xts_0.11-2                prettyunits_1.1.1        
#>  [23] jpeg_0.1-8.1              colorspace_1.4-1         
#>  [25] rvest_0.3.5               haven_2.2.0              
#>  [27] xfun_0.12                 callr_3.4.0              
#>  [29] crayon_1.3.4              jsonlite_1.6.1           
#>  [31] lme4_1.1-21               survival_3.1-8           
#>  [33] zoo_1.8-6                 glue_1.3.1               
#>  [35] gtable_0.3.0              webshot_0.5.2            
#>  [37] pkgbuild_1.0.6            rstan_2.19.2             
#>  [39] abind_1.4-5               scales_1.1.0             
#>  [41] DBI_1.1.0                 miniUI_0.1.1.1           
#>  [43] htmlTable_1.13.3          viridisLite_0.3.0        
#>  [45] xtable_1.8-4              progress_1.2.2           
#>  [47] HDInterval_0.2.0          ggstance_0.3.3           
#>  [49] foreign_0.8-72            Formula_1.2-3            
#>  [51] stats4_3.6.2              StanHeaders_2.19.0       
#>  [53] DT_0.11                   htmlwidgets_1.5.1        
#>  [55] httr_1.4.1                threejs_0.3.1            
#>  [57] arrayhelpers_1.0-20160527 RColorBrewer_1.1-2       
#>  [59] ellipsis_0.3.0            acepack_1.4.1            
#>  [61] pkgconfig_2.0.3           reshape_0.8.8            
#>  [63] loo_2.2.0                 farver_2.0.3             
#>  [65] nnet_7.3-12               dbplyr_1.4.2             
#>  [67] labeling_0.3              tidyselect_1.0.0         
#>  [69] rlang_0.4.5               reshape2_1.4.3           
#>  [71] later_1.0.0               munsell_0.5.0            
#>  [73] cellranger_1.1.0          tools_3.6.2              
#>  [75] cli_2.0.0                 generics_0.0.2           
#>  [77] gifski_0.8.6              broom_0.5.3              
#>  [79] ggridges_0.5.1            evaluate_0.14            
#>  [81] fastmap_1.0.1             yaml_2.2.1               
#>  [83] processx_3.4.2            fs_1.3.1                 
#>  [85] nlme_3.1-142              mime_0.9                 
#>  [87] mvnfast_0.2.5             rstanarm_2.19.2          
#>  [89] xml2_1.2.2                compiler_3.6.2           
#>  [91] shinythemes_1.1.2         rstudioapi_0.11          
#>  [93] png_0.1-7                 reprex_0.3.0             
#>  [95] tweenr_1.0.1              stringi_1.4.6            
#>  [97] highr_0.8                 ps_1.3.2                 
#>  [99] Brobdingnag_1.2-6         lattice_0.20-38          
#> [101] Matrix_1.2-18             nloptr_1.2.1             
#> [103] markdown_1.1              shinyjs_1.0              
#> [105] vctrs_0.2.3               pillar_1.4.3             
#> [107] lifecycle_0.1.0           bridgesampling_0.7-2     
#> [109] data.table_1.12.8         httpuv_1.5.2             
#> [111] R6_2.4.1                  latticeExtra_0.6-29      
#> [113] bookdown_0.16             promises_1.1.0           
#> [115] gridExtra_2.3             codetools_0.2-16         
#> [117] boot_1.3-23               MASS_7.3-51.4            
#> [119] colourpicker_1.0          gtools_3.8.1             
#> [121] assertthat_0.2.1          withr_2.1.2              
#> [123] shinystan_2.5.0           parallel_3.6.2           
#> [125] hms_0.5.3                 rpart_4.1-15             
#> [127] grid_3.6.2                minqa_1.2.4              
#> [129] coda_0.19-3               snakecase_0.11.0         
#> [131] rmarkdown_2.0             shiny_1.4.0              
#> [133] lubridate_1.7.4           base64enc_0.1-3          
#> [135] dygraphs_1.1.1.6